Statistical Properties of a 2D Granular Material Subjected to Cyclic Shear 
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This work focuses on the evolution of structure and stress for an experimental system of 2D photoelastic 
particles that is subjected to multiple cycles of pure shear Throughout this process, we determine the contact 
network and the contact forces using particle tracking and photoelastic techniques. These data yield the fabric 
and stress tensors and the distributions of contact forces in the normal and tangential directions. We then find 
that there is, to a reasonable approximation, a functional relation between the system pressure, P, and the mean 
contact number, Z. This relationship applies to the shear stress t, except for the strains in the immediate vicinity 
of the contact network reversal. By contrast, quantities such as P, t and Z are strongly hysteretic functions of 
the strain, £. We find that the distributions of normal and tangential forces, when expressed in terms of the 
appropriate means, are essentially independent of strain. We close by analyzing a subset of shear data in terms 
of strong and weak force networks. 
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INTRODUCTION 

In this work, we describe experiments that probe the mi- 
croscopic properties of sheared granular materials, with an 
eye towards understanding the statistical properties and small- 
scale phenomena which strongly influence larger scale behav- 
ior. The application of shear leads to the evolution of a strong 
force network, as shown in Fig. [T] sometimes referred to as 
force chains. These mesoscopic structures are filamentary net- 
works that carry forces at or above the mean, and that extend, 
in the case of shear, over distances of a few to perhaps many 
tens of grains. During shear, the force network evolves, with 
force chains strengthening, and then ultimately failing. 

The present studies explore the structural evolution of sys- 
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tern during shear by means of fabric, stress and related tensors. 
Associated with the evolution of the fabric and stress tensors 
are a number of complex phenomena, including shear bands, 
particle rotation, failure and buckling of the force chains, 
among other effects. We have recently shown the importance 
of rotation for the failure of force chains, particularly in shear 
bands 1 1]. In general, all quantities measured here show fluc- 
tuations, and, of course sensitivity to the direction of the shear 
(forward and reverse). In particular, when the shear direction 
is switched, the system undergoes structural reconstruction, 
causing changes in the average contact number, the mean ori- 
entation of the contacts, and the stress tensor. 



We focus on a path corresponding to pure shear strain, 
starting from a packing fraction where there is no observ- 
able stress. As we strain the system, the detected stresses 
and mean contact number Z increase, and the system reaches 
a jammed state for Z's above Z ~ 3. As we further deform 
the system, including reversal of the strain, Z tends to re- 
main at or above 3 for much of the time. Throughout, mul- 
tiple shear cycles, the packing fraction (|) remains at a fixed 
value, (|) = 0.758, that is below the observed jamming value 
for isotropic compression!^. When the shear strain is re- 
versed, the original force network largely vanishes, and a new 
strong network forms. This process is strongly hysteretic in 
the strain, but we find that the stresses can be characterized 
rather well in terms of the system-averaged contact number, 
Z. 



FIG. 1: Sequence of photoelastic images showing the evolution of 
the force chains as the system is sheared in the forward (images a 
and b), then the reverse direction (image c and d). These four images 
are chosen at different steps from the 1st shear cycle, with £ = 0.033, 
0.267, 0.267, and 0.033, respectively. The axial strain e is defined 
below. In images b and c, the sidewall of the biax has moved into 
view, creating dark bands at the bottom of the images. 



In the remainder of this work, we first describe basic fea- 
tures of the experimental techniques. We then present results 
from cyclic pure shear experiments, in which we explore the 
structural and stress changes within each shear cycle and in 
particular during shear reversals. We then analyze the force 
network in terms of strong and weak components. 
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EXPERIMENTAL TECHNIQUES 

The experiments described here use a 'biaxial' apparatus, 
as depicted in the upper parts of Fig.|2l This device allows us 
to deform a rectangular sample of particles into any other de- 
sired rectangular shape, hence apply pure shear (compression 
in one direction, but equal dilation in the other), uniaxial com- 
pression, isotropic compression, etc. Here, we focus uniquely 
on shear deformations, which maintain a fixed area for the 
system. 

The studies are carried out effectively in 2D by using disks 
which are made of photoelastic material. When under stress, 
and when viewed between crossed polarizers, photoelastic 
materials exhibit a series of light and dark bands, as in the 
bottom left image of Fig. |2l These bands encode the detailed 
stress within each particle, and these stresses are in turn, de- 
termined by the forces at contacts on each particle. In the 
past, this technique has been used in several different stud- 
ies yl |4[]. What makes our current approach unique is that, 
for large collections of particles, we solve the inverse prob- 
lem which starts from the photoelastic image and yields as 
output, the inter-particle contact forces. More details have 
been given elsewhere HI HQ]. The particles are also typically 
marked with a small bar, which allows us to track the rota- 
tion and displacement of individual particles. In our current 
tracking approach, the bars are drawn on with fluorescent ink, 
which is invisible under ordinary light but glows strongly un- 
der UV light. In this way, it is possible to have both photoelas- 
tic images for force/stress measurement and separate images 
for tracking rotation and displacement of a given set of par- 
ticles, without mutual interference. In earlier versions of this 
approach , we used soUd black bars drawn on 

the particles. In this case, we imaged separate sets of particles 
for determining forces and for tracking motion. 

The initial boundaries of the system form a square filled 
with 1568 bi-disperse photoelastic disks at a packing fraction 
(|) = 0.758. There are roughly 80% smaller particles having 
a diameter of 0.74 cm, and 20% larger particles having a di- 
ameter of 0.86 cm. The initial state is prepared as close to 
isotropic as possible and is stress-free. The system is then 
subjected to shear by compression along the y-direction and 
expansion along the x-direction, as in Fig.|2](top-left), keeping 
the system area constant. Since the area is fixed, the deforma- 
tion can be simply defined using the strain £ along the x-axis 
with e = {x — X(j) / XQ. Here, xq is the initial size of the square. 
Once a maximum deformation t,„ax is reached, shear is re- 
versed by compression along the jc-axis and expansion along 
the y-axis. After this first shear reversal, the shear continues 
until the system domain returns to a square, and then deforms 
in the negative strain direction. Once e reaches a minimum 
Emm < 0, a second shear reversal is applied, eventually return- 
ing the system domain to a square with e = 0. This completes 
one shear cycle. The second shear cycle continues from the 
final state of the first shear cycle and the same procedure is 
applied for a total of six shear cycles. Note that the actual val- 



TABLE I: A list of e,„„v and for different shear cycles. 
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FIG. 2: Top-left: Sketch of top view of the experimental apparatus, a 
2D 'biax', consisting of pairs of facing boundaries that can be moved 
precisely under computer control so as to produce desired strains. 
Particles rest on a smooth slippery sheet of Plexiglas and are confined 
laterally by the walls of the biax. Strains are applied quasi-statically, 
in small discrete steps. Top-right: Side view of apparatus. Imaging 
is carried out by a camera mounted above the biax, and for each step, 
we obtain three images: one with crossed polarizers (bottom left), 
one without polarizers (bottom center), and one without polarizers 
but with UV illumination (bottom right). 

ues of e,nax and e„„„ are different for each shear cycle and the 
possible extreme values for these are determined by the spatial 
limit of the apparatus. A list of £,„,„ and £,nax are summarized 
in Table |T] The whole shear process is carried out in small 
incremental quasi-static steps. From one step to the next, e 
increases or decreases by a small amount 5e = ±3.3 x 10^^, 
depending on the shear direction. After each step, the motion 
is paused and images are acquired. The three images in the 
bottom row Fig.|2]show close-ups of the three different image 
types. The left-most of these is taken with polarizers in place, 
the middle is without polarizers and with ordinary light, and 
the right-most is without polarizers and with UV light. 

Before we turn to detailed results we note an experimen- 
tal issue of importance. During the parts of the cycle where 
the overall stresses in the system are low, the photoelastic re- 
sponse at some contacts falls below our limit of resolution. 
Because, as discussed below, the experiments indicate a dis- 
tribution of normal contact forces of the form 

P{F„) = {F„)-\f {¥„/{¥„)), (1) 

where / is to a reasonable approximation, the same function 
for all mean forces, we can estimate the number of the missed 
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FIG. 3: Sketch illustrating the notation for calculating the fabric, 
force-moment and stress tensors. 

contacts reasonably well. Non-zero contact forces below our 
experimental resolution also affect our measurements of stress 
and, in particular, P. However, the effect on stress components 
is much lower, since the contact forces appear linearly in the 
appropriate sums. 

We expect that we miss a fraction 

r^P{Fn)dF„/ [ P{F„)dF„ (2) 
Jo Jo 

of contacts, where F^ is a small known cut-off force, roughly 
the weight of a particle, below which we cannot detect the 
photoelastic response. This means that the measured Z's are 
lower than their true values by 

/ P{F„)dF„/ I P{Fn)dF„. (3) 
Jfc Jo 

We also underestimate the pressure by a factor of 

/ F„P{F„)dF„/ f F„P{F„)dF„. (4) 
Jiv Jo 

In this last expression, we assume that all particles have the 
same radius, which is a reasonably good assumption. In this 
regard, the correction can be as much as 15% in Z very near 
jamming, but then becomes negligible for Z a bit above 3.0. 
The coiTection to P is much smaller, only 1-2% close to the 
jamming transition. In order to simplify the correction, we 
assume that the force distribution is an exponential. This is 
roughly right, and produces a reasonable correction, given the 
statistical variability of the data. The results of P and Z pre- 
sented in this paper have been corrected accordingly. We use 
the same correction factor for shear stresses, x, although in 
this case, the correction is not a rigorous. 

EXPERIMENTAL RESULTS 

We are concerned with the evolution of the force and con- 
tact networks. Both are typically strongly anisotropic, and the 
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FIG. 4: Evolution of the average contact number Z vs. step number 
Ns- Each shear cycle is colored differently. Arrows indicate the steps 
where the shear direction is switched. 
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FIG. 5: Mean contact number Z vs. strain e for the second cycle. 

direction of the anisotropy switches quickly when the direc- 
tion of strain is reversed. The force anisotropy is evident in 
Fig. [T] which shows representative photoelastic images dur- 
ing different phases of a single cycle. The structural changes 
of the contact network during cyclic shear are also strongly 
anisotropic. These can be captured by the fabric tensor, Rij, 
defined as 

<:=lt = l 

Here, the summation and include only non-rattler disks, and 
as illustrated in Fig. [3] c/^ is the number of contacts on disk k, 
and nli^ is the Zth component of the unit branch vector point- 
ing from the center of the disk to a contact c. We consider 
a rattler disk to have less than two detectable contacts. The 
average contact number Z is simply the trace of the fabric ten- 
sor Rjj. The principal eigendirection of Rij is also a useful 
measure of the prevailing orientation of the force network. 

Figs. |4] and |5] show how Z changes in each shear cycle as 
a function of step number Ns for all cycles, and as a function 
of strain £ for one cycle. In Fig. 21 the individual shear cy- 
cles are distinguished using different colors. We will maintain 
this color scheme throughout to identify the various cycles. 
Arrows in Fig.|4]indicate the shear reversals. Z fluctuates be- 
tween a minimum value of around 2.5 and a maximum value 
of around 3.5. The first shear cycle, the red curve, begins 
with a nearly stress-free and isotropic state. The force chains 
build up steadily. As a consequence, Z increases as more force 
chains develop. Z barely exceeds 3 before the first shear re- 



4 



1 

0.8 
0.6 

L 

0.4 
0.2 

Q 



200 400 600 800 1000 1200 1400 

Step 



200 
150 
<?J 1001^' 
50 

















r 
















^1 



200 400 600 800 1000 1200 1400 

Step 



FIG. 6: Evolution of ttie ratio, p, vs. step number, Ns. p is defined as 
tlie number of non-rattler particles over the total number of particles. 



versal. Immediately after the reversal, Z drops. The relatively 
rapid decrease of Z after a switch of the shear direction is 
common to all shear reversals because the force network/force 
chains switch direction during this transition. With continued 
strain after a reversal, Z again increases as a new strong net- 
work, orthogonal to its predecessor, emerges. Note that Z = 3 
is the nominal isostatic point for frictional particles in two di- 
mensions. Fig. m shows that the system remains mostly in 
a jammed state after the first cycle, but that it may leave a 
jammed state briefly after shear reversal. As we show below, 
shear bands form in response to the shear. Hence, these states 
are not spatially homogeneous. However, it is the presence of 
a mechanically rigid but dilated region in the shear bands that 
allows for jammed states at (|)'s lower than the isotropic value. 
We also emphasize that Z is strongly hysteretic when viewed 
with respect to strain. To demonstrate this point, we show one 
shear cycle, the second, as a function of strain in Fig.|5] 

At various steps, the fraction, p, of detectable non-rattler 
particles also fluctuates, as displayed in Fig. |6] This ratio 
changes from 0.2 up to 0.90. Both Fig. |4] and Fig. |6] show 
similar trends, although p is noisier at the step where shear is 
reversed. 

We characterize the mean anisotropy of the contact network 
in terms of i3^, the system-averaged value of the angle between 
the eigenvector of the maximum eigenvalue of R and the x 
axis. Here, we restrict < < 180". Figure |7] shows that 
'& switches quickly, shortly after each strain reversal. That 
is, after a very small strain, aligns with the compressive 
direction. In order to see how quickly the angle changes after 
a shear reversal, we have plotted on a much finer scale. The 
results are presented in Fig. [8] where the graphs are organized 
from top to bottom, as a function of step number, N^- The 
typical number of steps required for the readjustment of the 
orientation varies from ANs ~ 10. However, we note that there 
is often a lag between the strain reversal, and the switch in "d, 
which may occur in only a few strain steps. 

The stress tensor a,j and the force moment tensor, Oij pro- 
vide additional measures of anisotropy, in this case for the 
forces. We define a local force moment tensor as 



(6) 



FIG. 7: Fabric orientation angle ■& vs. step number, Ns- t3- is de- 
fined as the absolute angle between the eigenvector of the maximum 
eigenvalue of the fabric tensor R and the x axis. This angle measures 
the dominant contact orientation 
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FIG. 8: Fabric orientation angle i5 vs. step number Ns on a fine scale 
near shear reversals. Arrows in each graph indicate the beginning of 
the shear reversals. The readjustment of t3- after a shear reversal takes 
about 10 to 20 steps. 
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FIG. 9: Pressure P vs. step number Ns. 
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FIG. 10: Shear stress x vs. step number A'j. 
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FIG. 1 1: 8P (top) and 8t (bottom) vs. step number Ns- 8P and 5t are 
respective differences between two neighboring steps of P and t as 
given m Fig.|9]and Fig.[lO] 



FIG. 14: Pressure, P vs. strain, e, for the second cycle, showing 
strong hysteresis. 
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FIG. 12: PDFs of 8p„ (circles), and 5t„ (squares). 8p„ and 5t„ are 
made dimensionless through normalization by their corresponding 
standard deviations. 



The globally averaged stress tensor is then 



N 
k=l 



(7) 



Here, A is the system area; A^, Ck, i, r^^ and j have the same 
meaning as in the expression of Rjj (e.g. Fig. O. f-j^ is the 
ith component of the contact force on particle k at contact c. 
The two eigenvalues of the stress tensor are Oi and 02, where 
ai < 02 by definition. The pressure is then P — j{oi +O2) 
and the shear stress is x = 5 {02 — c?i )• 
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FIG. 13: Power spectra, on log-log scales, for 8p„, (inset) and 5t„ 
computed from the data given in Fig. [9] Fig. [TO] vs. frequency cor- 
responding to inverse step number, where a frequency of 1.0 corre- 
sponds to a 1/1000 steps. The solid line is a guide to the eye, and 
corresponds to a power law with an exponent of —2. 



FIG. 15: Shear stress x vs. strain, e for the second cycle. 



Figure|9]and Fig.[TO]show the evolution of the pressure and 
the shear stress versus step number, Ns, and Figs. [T4l and [TS] 
show P and T over the second cycle vs. strain. Both P and 
T vary significantly with strain, but unlike iJ, both quantities 
evolve steadily, modulo some significant fluctuations, up to 
the maximum of a given cycle. The fluctuations evident in the 
stress evolution are due to failure events, which are associated 
with a continual collapse of old and formation of new force 
chains. To show the scale of the fluctuations more clearly, 
we present in Fig. [TT] the changes bP and 5t in P and T, re- 
spectively, between successive steps. The curves exhibit ran- 
dom spikes in the positive and negative directions. It is per- 
haps worth emphasizing that the fluctuations, even for these 
system-averaged stress differences, can be large relative to 
the locally averaged step sizes, which cannot be distinguished 
from zero on these plots. Figure[T2]shows the probability dis- 
tribution functions (PDF's) of 5P and 5t, after normalizing by 
their respective standard deviations. The two PDF's for dp„ 
and 5t„ are virtually indistinguishable, and decay exponen- 
tially for positive and negative values, with a somewhat more 
rapid decay for positive values. Sharp drops of P and T occur 
at steps where the shear is reversed. Each drop is associated 
with the release of stored mechanical energy. Interesting is- 
sues include the nature of the energy dissipated, and the rela- 
tion to the evolution of the force network (|1, 



m\lM)- For 



instance, recently, we have found that buckling of the force 
chains yj] is an important mechanism leading to the loss of 
energy in shear. We will address these questions in future 
work. 
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FIG. 16: Pressure, P, vs. average contact number, Z. Different colors 
correspond to data from different shear cycles. Data points from 
the first shear cycle, shown in red, deviate slightly from other shear 
cycles. 



FIG. 19: Data for the distribution of normal forces, F„, expressed as 
Fn/{Fn) of all six shear cycles. We show data from all strain steps 
collectively for each cycle. Here, data for F„/{F„) are normalized by 
{F„) for the given step. 




Z 

FIG. 17: Shear stress, x vs. average contact number, Z. 

To further examine the fluctuations of P and T, we analyzed 
their power spectra from the curves in Figs.|9]and[T0l Here, we 
use step number, A'^- as a time-like variable, and (A(v/1000)^' 
as a frequency-like variable. The spectra vs. frequency vari- 
able, shown in Fig.[T3]are similar for P and T, and are typically 
broad-band. They suggest a power law decay with an expo- 
nent close to —2 for both P and T. Similar behavior for the 
high-frequency part of the spectrum has been reported in pre- 
vious experiments on continuously sheared 3D granular sys- 
tems il4ll . 

It is clear from Figs. |5] [141 and[T5]that the stresses and Z are 
hysteretic in the strain, i.e., that the strain does not provide a 
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FIG. 18: Shear stress, x vs. average contact number, Z after removal 
of data points from the first shear cycle and data points where the 
strong network direction, as measured by d, is switching directions, 
as in Fig. [8] 



unique characterization of a state. It is then interesting to ask 
whether there is some other quantity that better characterizes 
the nature of a given state, and in particular whether there is a 
relation among Z and P and T. For jamming of spherical (cir- 
cular in 2D) particles under isotropic stress conditions, the key 
control parameter is the density/packing fraction. Then, both 
Z and P are functions of (|), and consequently, P is a function 
of Z. Here, however, the packing fraction is constant. Yet, 
starting from an unjammed state, we arrive at a state which is 
jammed when we apply sufficient shear. In order to address 
what might control jamming in this case, we note that Figs.|5] 
[T4l and [15] show similar shapes. Drawing on the isotropic 
stress case, we ask whether a relation exists between Z and 
the stresses. 

Indeed, Fig. [16] shows that data for P vs. Z fall on a nearly 
common curve. In a similar fashion, x vs. Z falls on a nearly 
common curve. For T vs. Z, the relative scatter is higher 
However, there is a systematic part of the T vs. Z data that 
fall below the weight of the curve. These data correspond to 
relatively small ranges of strain following a reversal. As i5 
switches direction, the system passes through a more nearly 
isotropic state. In addition, T vs. Z data from the first cycle 
start from an isotropic state for which T = 0. For the first part 
of that cycle, the system retains some memory of its initial 
state. If the data from the first cycle and immediately after re- 
versals are removed, the results for x vs. Z yield a collapse 
that is comparable to that for P vs. Z as seen in Fig. [18] 
Although P and Z do collapse, the spread of the data points 
around the curve is still quite big. This spread reflects statis- 
tical fluctuations, which as shown above, can be large even 
for the system-averaged P and T. An interesting observation 
is that except for the switching regimes, in the mean, the ra- 
tio z/P = constant. In this case, the relation follows because 
T and P are separately (essentially) linear functions of Z and 
both vanish at the common value Z ~ 3 where the system first 
jams. 

Additional statistical measures include the distribution of 
contact forces, including the normal force distribution given 
in Fig.[T9]and the tangential force distribution given in Fig. [20] 
In these figures, we have organized the data by cycle number. 
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FIG. 20: Data for the distribution of normalized tangential forces, 
Ft/{F„) of all six shear cycles. As for the distribution of normal 
forces, the statistics are combined for all steps within each cycle, 
where the normalization, {F„) is made for each step. 
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FIG. 21: Data for the distribution of normal forces, F„, expressed 
as F„/{F„) in the first shear cycle, (a) distribution P(F„/{F„)) for 
forward shear, e > 0. (b) P{Fn/ {Fn)) for reverse shear when e > 0. 
(c) P{F„/{F„)) for reverse shear when e < 0. (d) P{F„/{F„)) for 
forward shear when £ < 0. To improve statistics, each data point on 
the graph includes a set of F„/ {F„) from 10 neighboring steps in most 
cases. A summary of each data point and its corresponding steps and 
strains can be found in Table Ull 



FIG. 22: Data for the distribution of tangential forces. Ft, expressed 
as Ff/{F„) in the first shear cycle, (a) distribution P{F,/{F„)) of the 
forward shear when the strain e > 0. (b) P{F,/{F„)) of the inverse 
shear when £ > 0. (c) P{F, / (F„) ) of the inverse shear when e < 0. (d) 
P(F,/{Fn)) of the forward shear when £ < 0. To improve statistics, 
each data point on the graph includes a set of Ft/{F„) from ~10 
neighboring steps. See details in Table Ull 



differ somewhat, but this is to be expected because the statis- 
tics are Hmited there. P{Ft/{Fn)) also shows an exponen- 
tial decay, except that the tangential force distribution, as a 
function of P{Ft/ {F„)), decays faster than that for the normal 
forces. However, this is simply due to the choice of normal- 
ization. That is, if the tangential forces were normalized by 
(Fi), then the rate of exponential fall-off would be comparable 
to that for the normal force distributions. For small Fn/ {Fn) 
and Ft/{F„), both distributions fall below exponentials. Al- 
though some of the fall-off is due to the experimental lower 
limit of force detection, we believe that this is a relatively mi- 
nor effect. Specifically, the shape of the distribution is not 
particularly sensitive to the mean force. 



and we have combined data for different steps within a cycle. 
Here, we normalized the data for F„ or Ft for each step by the 
the mean normal force {Fn) at that step. To justify that this 
is legitimate, we plot the distributions of Fn/ {F„) and Ft/ {Fn) 
for different strain steps within the first shear cycle in Figs. 1211 
and |22] Several neighboring steps, ten for most data points 
and three to five for data near a reversal, are combined for each 
data point on the plots. The details are summarized in Table|II] 
Some strain steps may have a longer tail than others but their 
general shapes are more or less similar, and in particular, there 
is no systematic difference between distributions for different 
steps. As always, the tails show bigger scatter due to limited 
statistics. 

The distributions of normal forces show a common form 
consisting of a nearly exponential fall-off at large Fn/ {F„) and 
a peak at low Fn/{F„). At the extreme tails, the distributions 



CONNECTION TO SHEAR LOCALIZATION AND FORCE 
CHAIN EVOLUTION 

We next explore possible connections between the hystere- 
sis observed in these experiments and two defining aspects of 
material behavior under shear: shear banding and stick-slip 
(e.g. [il2i Il5[l ). These two mechanisms are related and are 
both governed by force chain/force network evolution. Note 
that we use the term 'stick-slip' to signify fluctuations in the 
macroscopic stress, in particular, that of the stress ratio, z/P. 
Specifically, stick and slip events are periods where the stress 
ratio increases and decreases, respectively, with increasing 
strain. A slip event is most often due to the collapse of force 
chains by buckling and hence is accompanied by the release 
of stored energy, accumulated in force chains during the pre- 
ceding stick event. 

In what follows, we examine the evolution of the force 
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TABLE II: Summary of step numbers and strains for each group la 
beled data,- in Figs. l21|[22l 



Panel 


data; 


steps 


£ 


a 


1 


2-11 


0.0067 ~ 0.0367 




2 


12-21 


0.04 ~ 0.07 




3 


22-31 


0.0733 ~ 0.1033 




4 


32-41 


0.1067 ~ 0.1367 




5 


42-51 


0.14 ~ 0.17 




6 


52-63 


0.1733 ~ 0.21 




7 


64-73 


0.2133 ~ 0.2433 




8 


74-83 


0.2467 ~ 0.2767 




9 


84-86 


0.28 ~ 0.2867 


b 


9 


88-97 


0.28 ~ 0.25 




8 


98-107 


0.2467 ~ 0.2167 




7 


108-117 


0.2133 ~ 0.1833 




6 


118-127 


0.18 ~ 0.15 




5 


128-137 


0.1467 ~ 0.1167 




4 


138-147 


0.1133 ~ 0.0833 




3 


148-157 


0.08 ~ 0.05 




2 


158-167 


0.0467 ~ 0.0167 




1 


168-177 


0.0133 ~ -0.0167 


c 


1 


178-187 


-0.02 ~ -0.05 




2 


188-197 


-0.0533 ~ -0.0833 




3 


198-207 


-0.0867 ~ -0.1167 




4 


208-217 


-0.12 ~ -0.15 


d 


5 


218-227 


-0.1467 ~ -0.1167 




4 


228-237 


-0.1133 ~ -0.0833 




3 


238-247 


-0.08 ~ -0.05 




2 


248-257 


-0.0467 ~ -0.0167 




1 


258-302 


-0.0133 ~0 



network and its effect on the macroscopic stress. We con- 
fine our attention to the different time steps for the forward 
shear. These data are from a run with a slightly higher den- 
sity, (|) = 0.76. As shown earlier, the force network evolution 
in the forward and reverse shear cycles exhibit fairly univer- 
sal statistics. To proceed, we employ two algorithms. The 
first identifies force chain particles via a so-called particle 
load vector. For each particle, the local force moment tensor, 
Ojj, as defined earlier, is computed. The largest eigenvalue 
of this tensor and its associated eigenvector are then used to 
define, respectively, the magnitude and direction of the parti- 
cle load vector. The direction of force transmission is dictated 
by the direction of the particle load vector. Groups of parti- 
cles whose particle load vectors line up within a predefined 
narrow angular range, and whose particle load magnitude is 
above the global average value, constitute a force chain. This 
procedure has been incorporated into an algorithm that takes 
contact force data as the known input, and provides the force 
chain particles, and hence the force chain particle network, as 
the output. Complete details of this procedure and its associ- 
ated algorithm are provided elsewhere (idi 17 ] . 

The second algorithm was developed for the purposes of 
identifying parts of the force chain particle network that have 
undergone buckling, i.e. buckled force chain "BFC" segments 




FIG. 23: Displacement fields showing shear bands in forward and 
reverse shear. These four fields correspond to the small particle dis- 
placement right after applying small deformations to the four images 
shown in Fig.[T| The red lines drawn in (b) (d) are a guide to the eye, 
indicating the regions of shear localization. Note that in parts (c) and 
(d) one or more side walls occupies part of the image. 




1 1211 . A strain interval of interest, [e'^,e^] is chosen: for exam- 



ple, that which spans a drop in stress ratio or an "unjamming 
event", or a single time step in a DEM simulation. A set of 



FIG. 24: Force chains in the shear band regime. This image is a slice 
of the main diagonal regime of image (b) in Fig.[T] Here the image 
is rotated to make the original diagonal line horizontal. 



three filters is then applied: (a) eliminate all particles not in 
force chains at e'^; (b) out of those remaining, eliminate those 
which have not decreased in potential energy; (c) out of those 
remaining, identify and isolate all 3-particle segments which 
have buckled. To determine if a segment has buckled, we con- 
sider the angle between the branch vectors from the central 
particle to the two outer particles. The decrease in this an- 
gle over the interval in question is defined as being twice the 
buckling angle, 9^,. Then, a buckling segment is simply one 
where 9;, > 0. The set of particles remaining after all three fil- 
ters have been applied is the set refeiTed to hereafter as BFCs. 
Later, we consider populations of BFCs in distinct subsets, 
where each subset is distinguished by a predefined nonzero 
buckling threshold 9^* that member BFCs must satisfy over the 
given strain interval. Complete details of this entire procedure 
and associated algorithm are provided elsewhere lfl2ll . 

The specimen deforms in the presence of a shear band, as 
shown in Fig. |23] The band is backward inclined in the for- 
ward shear and forward inclined in the reverse shear. As the 
material is sheared in a given direction, two triangular blocks 
slide over the shear band in opposing directions, effectively 
leading to simple shear across the shear band BlSll . 

Although the material is dilated in the shear band, the force 
chains pass right through these bands, without substantial 
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FIG. 25: Evolution with axial strain, £22, of the population of force 
chain particles, A'^''. 




FIG. 26: Spatial distribution of force chain particles undergoing 
buckling (red) across the strain interval £22 = 0.060 to 0.077, to- 
gether with their confining neighbors (blue). A buckling threshold of 
Q* = 2° is used. Data shown here correspond to the central portion 
of the sample in the early stages of shear band development. We note 
that the width of the shear band narrows for higher strains. 




FIG. 27: Evolution with axial strain, £22, of population of buckled 
force chain segments, N^^^ , for various buckling thresholds, 6^. 
Green, blue and red lines correspond to 9^ = 1°, 2° and 3°, respec- 
tively. Also shown is stress ratio (black). 




0.00 



0.00 



0.02 



0.04 



0.06 



0.08 



FIG. 28: Evolution with axial strain, £22, of the relative rotation at 
contacts per step, 0, for various types of contacts. Black, green, blue 
and red correspond to EC (force chain) to EC, EC to WN (weak net- 
work), WN to WN and BEC (buckling force chain) to BEG contacts, 
respectively. The unit of 9 is radians. 



changes, due to force balance (Fig. 1241). During 'stick events', 
regimes of strain where the force network is stable, the pri- 
mary force chains provide the major resistance to motion in 
the compressive direction. Weak secondary force chains still 
exist in the dilation direction and serve to 'prop up' the pri- 
mary force chains. Primary force chains, laterally confined by 
weak network neighbors, are subject to axial compression and 
often fail via buckling. Secondary force chains, being in the 
direction of shear, tend to fail by extension. 

As mentioned earlier (see discussion around Figs. |9] and 
[Tol l, there is evidence of local failure events in the force net- 
work throughout the loading, both in the forward and reverse 
shear. These failure events are due to the continual collapse of 
old and formation of new force chains. The failures are con- 
centrated mainly in the shear band where the mode of defor- 
mation of the material is essentially one of simple shear MIST . 
To unravel the mechanisms behind the hysteresis in the nor- 
mal and shear components of stress, we examine the contribu- 
tion to these stresses from force chains versus weak network 
particles. Recall that we confine our attention to the forward 
shear in the first shear cycle. In Fig. |25] we show the strain 
evolution of the population of force chain particles. As ex- 
pected, there is an initial increase in the population as load is 



increased, followed by a near constant value for large strains. 
Note here, that the threshold value for the particle load vector 
(i.e. the global average) of each particle classified as a force 
chain is increasing with strain. 



In light of Figs. 123b andl23tl. the relative rotations between 
force chains are particularly relevant here. Discrete element 
simulations and photoelastic disk experiments of Oda and co- 
workers (e.g. 13 3; 21. 22]) and by Veje et al. and Utter et 
al.|8, 9] have shown that particle rotations concentrate in the 
shear band. More recent studies confirm this, but addition- 
ally found that rotational motions dominate during slip events 
inside the shear band (|[ll[l2, 13]). These slip events are gov- 
erned by the failure of force chains by buckling. The location 
and deformation periods during which relatively large parti- 
cle rotations occur do indeed coincide with the location and 
incidences of force chain buckling: see Fig.|26]and Fig. l27l 
Rotation is a key mechanism in force chain buckling as is evi- 
dent in Fig.|28]which shows that the greatest relative rotations 
are sustained at contacts between particles in buckling force 
chains. 
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CONCLUSIONS 

In this work, we have explored the evolution of force and 
contact networks for cyclic shear of a dense granular mate- 
rial. Starting from an initially unjammed low-density state, 
jamming occurs for sufficiently large shear. This is, in fact, 
associated with the phenomena of Reynolds dilatancy. In par- 
ticular, during the first shear cycle, the system reaches Z = 3 
around step 70 corresponding to a shear strain of e = 23.3% 
. During much of each cycle, the system is above a jamming 
threshold for which the mean contact number is Z ~ 3. Al- 
though the stress components, Z, etc. are strongly hysteretic 
in strain, we empirically find that P is strongly correlated with 
the mean contact number when the system is jammed. The 
same is true for x vs. Z except immediately after a strain rever- 
sal, when T can become small. Again, excluding regions im- 
mediately following reversals, the mean values for P{Z) and x 
approach zero linearly, within our resolution, as Z — » 3 from 
above. In this case, the ratio x/P is, on average, constant. 

We have apphed a correction to the data to account for con- 
tacts which are below the experimental detection threshold. 
This correction was most important for very low stress states. 
For Z, the correction could reach ~ 15%, although it was also 
applied to the pressure and the shear stress, x. The correc- 
tion was implemented by assuming that the force distributions 
have (close to) universal exponential forms. A more accurate 
correction technique will be implemented in future work, but 
the present approximation is reasonable, given the statistical 
scatter of the data. We note that the jammed states are gen- 
erally inhomogeneous in the density/packing fraction, since 
these are characterized by shear bands where the material is 
locally dilated and where much of the motion occurs. 

We then analyzed the evolving force network in the for- 
ward shear using algorithms that distinguish contributions to 
macroscopic stress of particles from the strong and weak con- 
tact force network. Relatively large rotations develop during 
the buckling of force chains. These buckling events, which 
are present throughout the loading history, are primarily con- 
fined to the shear band and dominate during sUp events or 
periods where stress ratio decreases with increasing strain. A 
detailed analysis of such internal failure events by buckhng 
and their cormection to structural evolution of the force and 
contact networks in cycUc shear loading is the subject of an 



ongoing investigation. 
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